clear all 
set more off 
set maxvar 15000 
clear matrix


    use "$Mydirectory2/appendix_d/DM_sample_foranalysis.dta", clear

    levelsof byear, local(levels) 
    di "`levels'"
       
/*Notes: (1) The "common sample" = respondents with 
             non-missing actual income in both 
             generations and non-missing predicted 
             working father income.
*/

*---------------------------*
* Create ranked measures *
*---------------------------*
/*
  (1) Actual parental income, common sample
  (2) Predicted father income, common sample
*/
    preserve
        keep if comb_sample==1 
        bysort byear: egen N_byear_combsample = sum(comb_sample) //# of obs in the baseline sample in each birth year

        foreach var of varlist lfaminc1 lfaminc0 log_father_interpolated_retro { 
        
        if "`var'" == "lfaminc1" local rname0 "c_rank"
        if "`var'" == "lfaminc1" local ysel0 "ysel_c_rank"
        
        if "`var'" == "lfaminc0" local rname0 "MDp_rank"
        if "`var'" == "lfaminc0" local ysel0 "ysel_MDp_rank"

        if "`var'" == "log_father_interpolated_retro" local rname0 "EJp_rank"
        if "`var'" == "log_father_interpolated_retro" local ysel0 "ysel_EJp_rank"

        
            egen `rname0' = xtile(`var') if inrange(byear,1948,1964), by(byear) nq(100) weight(weight)
            replace `rname0' =. if N_byear_combsample<100
            
            qui:gen `ysel0'=.
            replace `ysel0'= byear if `rname0'!=.

        }

        assert c_rank==. if lfaminc1==. 
        assert MDp_rank==. if lfaminc0==. 
        assert EJp_rank==. if log_father_interpolated_retro==.
        
        label var c_rank "Rank adult child, family income, common sample"
        label var MDp_rank "Rank parent, D/M inc, common sample"
        label var EJp_rank "Rank parent, Jácome et al inc, common sample"


*-----------------------*
* Run regressions *
*-----------------------*
    foreach pvar in MDp_rank EJp_rank {

        if "`pvar'"=="MDp_rank" local name "MD_combsamp"
        if "`pvar'"=="EJp_rank" local name "EJ"

        matrix rank_`name'_EJreg = J(10,3,.)
        local i =0

        foreach y in `levels' {
            reg c_rank `pvar' if byear==`y' [w=weight], cluster(hhid)

            local i=`i'+1
            matrix rank_`name'_EJreg[`i',1] = `y'
            matrix rank_`name'_EJreg[`i',2] = _b[`pvar']
            matrix rank_`name'_EJreg[`i',3] = _se[`pvar']
        }
    }

    restore

*------------------------------------------------------*
*------------------------------------------------------*

*-----------------------------------------------*
* Make figures *
*-----------------------------------------------*

    //Figure 1: overlay (1) and (2)
    clear
    svmat rank_MD_combsamp_EJreg
    ren rank_MD_combsamp_EJreg1 year
    gen l95_MD_combsamp_EJreg = rank_MD_combsamp_EJreg2-1.96*rank_MD_combsamp_EJreg3
    gen u95_MD_combsamp_EJreg = rank_MD_combsamp_EJreg2+1.96*rank_MD_combsamp_EJreg3
    sort year
    tempfile rank_overlay_EJreg
    save `rank_overlay_EJreg'

    clear
    svmat rank_EJ_EJreg
    ren rank_EJ_EJreg1 year
    gen l95_EJ_EJreg = rank_EJ_EJreg2-1.96*rank_EJ_EJreg3
    gen u95_EJ_EJreg = rank_EJ_EJreg2+1.96*rank_EJ_EJreg3
    sort year

    merge 1:1 year using `rank_overlay_EJreg'
    drop _merge

    gen year2 = year + 0.15

    preserve 
        keep if year<=1960 //birth cohorts up to 1960
        #delimit ;
        twoway (connect rank_MD_combsamp_EJreg2 year, msymbol(circle_hollow) mcolor(blue) msize(medium) lcolor(blue) lwidth(0.25) lpat(solid)) (rcap u95_MD_combsamp_EJreg l95_MD_combsamp_EJreg year, lcolor(blue) lwidth(0.25))
               (connect rank_EJ_EJreg2 year2, msymbol(square) mcolor(stone*1.25) msize(medium) lcolor(stone*1.25) lwidth(0.25) lpat(solid)) (rcap u95_EJ_EJreg l95_EJ_EJreg year2, lcolor(stone*1.25) lwidth(0.25)),        
        graphregion(color(white)) legend(on ring(0) row(2) pos(11) order(1 "Actual parental income" 3 "Predicted parental income")) 
        xtitle(" " "Birth Year") ytitle("Rank-rank slope" " ", axis(1)) ylabel(0(0.2)0.6,nogrid axis(1)) xscale(range(1947.5 1953.5))
        ;
        #delimit cr
        graph export "$Mydirectory2/appendix_d/rank_overlay_mdej_EJreg_upto1960.pdf", as(pdf) replace
    restore
